Example - 2D Grid Interpolation

The xmsinterp library contains classes for performing 2d (x,y) linear, natural neighbor, and idw interpolation: InterpLinear and InterpIdw. Natural neighbor interpolation is performed using the linear interpolation class. 3d (x,y,z) interpolation can be performed using the idw class. Points (x,y) must be given to create the interpolation class.

We will learn how to use the interpolation classes by using a set of 2D points (x,y) with a measurement at each point. These points come from a csv file. Then we will generate a grid that bounds these points and use different interpolation options to visualize the data.

In [1]:
import param, panel as pn, holoviews as hv
from holoviews.streams import Pipe
from bokeh.document import Document
import numpy as np
import xmsinterp as interp
import pandas as pd
import xarray as xr
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import time
hv.extension('bokeh', 'matplotlib')
In [2]:
# make the plots bigger by default
plt.rcParams['figure.dpi'] = 120
# read in the data to do the interpolation
df = pd.read_csv("plumedat.csv")
plot = df.plot(kind='scatter', x='x', y='y', c='c', title='sampled concentration', cmap='rainbow')
# drop the id columns from the data frame so we can pass the remaining data as points to the interp class
df.drop(['id'], axis=1, inplace=True)

Generate grid for visualization

Get the min, max (x, y) of the data and generate a grid of points.

In [3]:
# get the (x,y) range of the data and expand by 10%
x_range = (df.loc[df['x'].idxmin()]['x'], df.loc[df['x'].idxmax()]['x'])
buf = (x_range[1] - x_range[0]) * .1
x_range = (x_range[0] - buf, x_range[1] + buf)
y_range = (df.loc[df['y'].idxmin()]['y'], df.loc[df['y'].idxmax()]['y'])
buf = (y_range[1] - y_range[0]) * .1
y_range = (y_range[0] - buf, y_range[1] + buf)

# adjust these values for a higher/lower resolution grid
number_of_x_points=600
number_of_y_points=400
x = np.arange(x_range[0], x_range[1]+1, (x_range[1]-x_range[0])/number_of_x_points)
y = np.arange(y_range[0], y_range[1]+1, (y_range[1]-y_range[0])/number_of_y_points)
grid = np.meshgrid(x, y)
pts = [x for x in zip(*(x.flat for x in grid))]
xx, yy = grid

Interpolate to the grid using linear interpolation and display the results with matplotlib.

In [5]:
linear = interp.interpolate.InterpLinear(df.values)
z = linear.interp_to_pts(pts)

zz = np.reshape(z, (len(y), len(x)))
# make the plots bigger by default
plt.rcParams['figure.dpi'] = 120
p = plt.contourf(xx,yy,zz)
cbar = plt.colorbar(p)

Here is an example visualizing with xarray.

In [6]:
xa = xr.DataArray(np.reshape(linear.interp_to_pts(pts), (len(y), len(x))))
xa.plot()
Out[6]:
<matplotlib.collections.QuadMesh at 0x7f0ff8bed390>

Interpolate using natural neighbor. The first example shows a "constant" nodal function and the second example shows a "quadratic" nodals function. A nodal function can better model trends that exist in the measured values.

In [7]:
linear.set_use_natural_neighbor()
z = linear.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz)
cbar = plt.colorbar(p)

Now do natural neighbor with a quadratic nodal function.

In [8]:
linear.set_use_natural_neighbor(nodal_func_type="quadratic", nd_func_pt_search_opt="nearest_pts")
z = linear.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz)
cbar = plt.colorbar(p)

Now we will use Inverse Distance Weighted (IDW) to interpolate to the grid. Notice that unlike Linear and Natural Neighbor interpolation, IDW interpolation can extrapolate beyond the bounds of the input data points.

In [10]:
idw = interp.interpolate.InterpIdw(df.values)
z = idw.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz)
cbar = plt.colorbar(p)

Now interpolate using a nodal function to try to better capture trends in the data set.

In [11]:
idw.set_nodal_function(nodal_func_type="gradient_plane")
z = idw.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz, 10)
cbar = plt.colorbar(p)

Notice, in the previous interpolation that we get negative values in some locations. If we assume that there should be no negatives values (like when measuring concentration) we can use the truncation options to adjust the interpolation.

In [14]:
idw.set_truncation_max_min(smax=150, smin=0)
z = idw.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz, 10)
cbar = plt.colorbar(p)

Interactive interpolation using param, holoviews, and panel

Below we define a param class for the interpolation options. Then we use panel to show widgets and a dynamic plot. The user can change the widget values and then push the "Compute new plot" button to update the plot using the currently specified interpolation options.

In [17]:
pipe = Pipe([])
bounds = x_range[0], y_range[0], x_range[1], y_range[1]

class InterpOptions(param.Parameterized):
    interpolation_option = param.ObjectSelector(default='idw', objects=['linear', 'natural_neighbor', 'idw'],
                                                precedence=1)
    nodal_function = param.ObjectSelector(default='constant', objects=['constant', 'gradient_plane', 'quadratic'],
                                          precedence=1.1)
    truncation = param.Boolean(default=True, precedence=2)
    truncation_maximum = param.Number(default=df.loc[df['c'].idxmax()]['c'], precedence=2.1)
    truncation_minimum = param.Number(default=df.loc[df['c'].idxmin()]['c'], precedence=2.2)
    compute_new_plot = param.Action(lambda x: x.update_hv_plot(), precedence=0)
    index = 0
    
    @param.depends('interpolation_option', watch=True)
    def _update_nodal_function(self):
        if self.interpolation_option == 'linear':
            self.params('nodal_function').precedence = -1
        else:
            self.params('nodal_function').precedence = 1.1

    @param.depends('truncation', watch=True)
    def _update_truncation_opts(self):
        if self.truncation:
            self.params('truncation_maximum').precedence = 2.1
            self.params('truncation_minimum').precedence = 2.2
        else:
            self.params('truncation_maximum').precedence = -1
            self.params('truncation_minimum').precedence = -1
    
    def generate_zz(self):
        # Get Interp Object
        interp_object = None
        if self.interpolation_option != 'idw':
            #interp_object = linear
            interp_object = interp.interpolate.InterpLinear(pts=df.values)
            if self.interpolation_option == 'natural_neighbor':
                self.index += 1
                interp_object.interp_to_pt((0,0)) # this will force triangle creation
                interp_object.set_use_natural_neighbor(nodal_func_type=self.nodal_function,
                                                       nd_func_pt_search_opt="nearest_pts")
        else:
            interp_object = interp.interpolate.InterpIdw(pts=df.values)
            interp_object.set_nodal_function(nodal_func_type=self.nodal_function)

        if self.truncation:
            interp_object.set_truncation_max_min(self.truncation_maximum, self.truncation_minimum)
        z = interp_object.interp_to_pts(pts)
        bob = np.reshape(z, (len(y), len(x)))
        return bob
        
    def update_hv_plot(self):
        pipe.send([self.generate_zz()])
        
    def get_image(self):
        return hv.Image(self.generate_zz(), bounds=bounds)
    
      
opts = InterpOptions(name="")
#opts.get_image()
In [18]:
%%output size=200
dmap = hv.DynamicMap(lambda **kw: opts.get_image(), streams=[pipe])
pn.Row(opts, dmap)
Out[18]: